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ABSTRACT 

We study the effect of radiative feedback on accretion onto intermediate mass black holes (IMBHs) using 
the hydrodynamical code ZEUS-MP with a radiative transfer algorithm. In this paper, the first of a series, 
we assume accretion from a uniformly dense gas with zero angular momentum and extremely low metallicity. 
Our ID and 2D simulations explore how X-ray and UV radiation emitted near the black hole regulates the gas 
supply from large scales. Both ID and 2D simulations show similar accretion rate and period between peaks 
in accretion, meaning that the hydro-instabilities that develop in 2D simulations do not affect the mean flow 
properties. We present a suite of simulations exploring accretion across a large parameter space, including 
different radiative efficiencies and radiation spectra, black hole masses, density and temperature, Too, of the 
neighboring gas. In agreement with previous studies we find regular oscillatory behavior of the accretion rate, 
with duty cycle ~ 6%, mean accretion rate 3% (Too/10 4 K) 2 5 of the Bondi rate and peak accretion ~ 10 
times the mean for ranging between 3000 K and 15000 K. We derive parametric formulas for the period 
between bursts, the mean accretion rate and the peak luminosity of the bursts and thus provide a formulation 
of how feedback regulated accretion operates. The temperature profile of the hot ionized gas is crucial in 
determining the accretion rate, while the period of the bursts is proportional to the mean size of the Stromgren 
sphere and we find qualitatively different modes of accretion in the high vs. low density regimes. We also find 
that softer spectrum of radiation produces higher mean accretion rate. However, it is still unclear what is the 
effect of a significant time delay between the accretion rate at our inner boundary and the output luminosity. 
Such a delay is expected in realistic cases with non-zero angular momentum and may affect the time-dependent 
phenomenology presented here. This study is a first step to model the growth of seed black holes in the early 
universe and to make a prediction of the number and the luminosity of ultra-luminous X-ray sources in galaxies 
produced by IMBHs accreting from the interstellar medium. 

Subject headings: accretion, accretion disks - black hole physics - hydrodynamics - radiative transfer - dark 
ages, reionization, first stars - methods: numerical 



1. INTRODUCTION 

The occurrence of gas accretion onto compact gravitating 
sources i s ubiquitous in the universe. The B ondi accretion 
formula (Bond i & Hoylel Il944t iBondil 119521) . despite the 
simplifying assumption of spherical symmetry, provides 
a fundamental tool for understanding the basic physics 
of the accretion process. Angular momentum of accreted 
gas, in nearly all realistic cases, leads to the formation 
of an accretion disk on scales comparable to or possibly 
much greater than the gravitational radius of the black hole, 
r g ~ GM/c 2 , thus breaking the assumption of spherical 
symmetry in the Bondi solution. However, the fueling of 
the disk from scales larger than the circularization radius 
r c ~ j 2 /GM, where j is the gas specific angular momen- 
tum, can be approximated by a quasi-radial inflow. Thus, 
assuming that numerical simulations resolve the sonic radius, 
r s , the resolved gas flow is quasi-spherical if r c <C r s . 
The Bondi formula, which links the accretion rate to the 
properties of the environment, such as the gas density and 
temperature, or Eddington-limited rate are often used in 
cosmological simulations to model the supply of gas to 



the a ccre tion disk from galactic scales (|Volonteri & Rees 
| 2005|: iDi Matteo. Colberg. Springel. Hernquist. & Siiacki 



2008; 



Pelupessv, Di Matteo, & Ciardi 



2007 



2008 



Greif, Johnson, Klessen, & Bromm 
Ulvarez. Wise. & Abelll2009l) . 

However, the Bondi formula is a crude estimation of the 
rate of gas supply to the accretion disk because it does not 
take into account the effect of accretion feedback loops 
on the surrounding environment. Radiation emitted by 
black holes origin ates from gra vitational potential energy 
of inflowing gas (Shapirol 119731) and a substantial amount 
of work has been performed to understand the simplest 
case of spherical accretion onto compact X-ray sources 
or quasars. Several authors have used hydrodynamical 
simulations to explore how feedback loops operates and 
whether they produce time-dependent or a steady accretion 
flows. A variety of feedback processes have been con- 
sidered: X-ray preheati ng, gas cooling, photo-heating and 
radiat ion press ure (lOstriker. Weaver 1 Y ahil. & McCra 
1971 ICowie. Ostriker & Starkl " ~Tl97" 
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Bisnovatyi-Kogan & Blinnikov 
1983[ iVitelld fl984 




1980; Krolik & London 



Wandel, Yahil, & Milgrom 
Milosavlievic. Bromm. Couch. & Oh 

120 lOt 



iNovak . Ostriker, & Ciotti 



2 



Park and Ricotti 



lOstriker. Choi. Ciotti. Novak. & Progal 120101) . Typically, 
the dominance of one process over the others depends on the 
black hole mass and the properties of the gas accreted by 
the black hole. The qualitative description of the problem 
is simple: gravitational potential energy is converted into 
other forms of energy such as UV and X-ray photons or 
jets, which act to reduce and reverse the gas inflow, either 
by heatin g the gas or by prod u cing momentum driven 
outflo ws dCiotti & Ostrikeri 120071: ICiotti. Ostriker & Progal 
I2009t iProgal 120071: iProga. Ostr iker. & Kurosawa! 120081) ? 
In general these feedback processes reduce the accre- 
tion r ate and thus the luminosity of the accr e ting black 
hole (lOstriker. Weaver. Yahil. & McCravl 119761: iBegelmanl 
119851: IRicotti. Ostriker & Mackl 120081) . Consequently, 
the time averaged accretion rate differs from Bondi's 
solution. There have been works on self-regulated 

accretion of supermassive black holes (SMBH s) at 

the center of elliptic a l galaxies (ISazonov et al.1 120051: 
ICiotti & Ostrikeri 120071 ICiotti et alj 120091) and radiation- 
driven axisymmetri c outflow in active galact i c nu- 
clei (IProgal 120071: IProga Ostriker. & Kurosawal 120081: 
Kurosawa, Proga, & Nagamine] 120091 : iKurosawa & Progal 
2009allbl) . However, far less has been done in order to quantify 
the simplest case of spherical accretion onto IMBHs as a 
functi on of the properties of the environment. Rece nt theoret- 
ical (|Milosavlievic, Bro mm. Couch. & Ohl 120094 hereafter 
MBCO09) and numerical dMilosavljevic. C ouch. & Bromml 
2009b, hereafter MCB09) works explore accretion of pro- 
togalactic gas onto IMBHs in the first galaxies. MCB09 
describes the accretion onto a 100 M Q black hole from 
protogalactic gas of density i%h,oc = 10 7 cm -3 and tem- 
perature Too = 10 4 K. Our study, which complements this 
recent numerical work, is a broader investigation of accretion 
onto IMBHs for a set of several simulations with a wide 
range of radiative efficiencies, black hole masses, densities 
and sound speeds of the ambient gas. Our aim is to use 
simulations to provide a physically motivated description of 
how radiation modifies the Bondi solution and provide an 
analytical formulation of the problem (see MBCO09). 

The results of the present study will help to better un- 
derstand the accretion luminositi es of IMBH s at high z 
and in the present-day universe dRicottil 120091) . Applica- 
tions of this work include studies on the origin of ultra- 
lumin ous X-ray sources (U LX s) dKrol ik. McKee, & Tarter 
19811 iKrolik & Kallmanl 119841: iRicottil 120071: 
Mack. Ostriker. & Ricottil 120071: IRicotti. Ostriker. & Mackl 
20081). the build up of an early X-ray background 
(IVenkatesan. Giroux. & Shulll 1200 It IRicotti & Ostriker 
20041: iMadau. Rees. Volonteri. Haardt. & Ohl 12004 



Ricotti, Ostriker, & Gnedin 2005) and growth of S MBHs 



IVolonteri. Haardt & Madaul 120031: IVolonteri & R ees 2005; 
Johns on & Bromml 120071: iPelupessy. Di Matteo. & Ciardil 
20071: lAlvarez. Wise. & Abell 120091) . For example, different 



scenarios h ave been propos ed for the formation of quasars 
at z ~ 6 (Fan et al. 2003): growth by mergers, accretion 
onto IMBHs, or direct formation of larger seed black 
holes from collapse of quasi-stars dCarr. Bond. & Arnef 



lett 

98; 



Frver, Wooslev. & Heger 2001; Begelman, Volonteri, & Ree; 
20061: IVolonteri. Lodato. & Nataraianl 120081 

Omukai, Schneider, & Haiman 2008; iRegan & Ha ehnelt 
20091: iMaver. Kazantzidis. Escala. & Callegaril 120101) that 
may form from metal free gas at the center of rare dark matter 



halos dOh & Ha iman 2002). Understanding the properties 
which determine the efficiency of self-regulated accretion 
onto IMBHs is important to estimate whether primordial 
black holes produced by Pop III stars c an accrete fast enough 
to become SMBHs by redsh ift z ~ 6 (IMadau & Reesll2001l: 
I Volonteri. Haardt. & Madaul 120031: lYoo & Miralda-Escudd 
20041: IVolonteri & Reesl 120051: Uohnson & Bromml 120071: 
Pelupessv. Di Matteo. & Ciardill2007t lAlvarez. Wise. & Abell 
2009). 

In this paper we focus on simulating accretion onto IMBH 
regulated by photo-heating feedback in ID and 2D hydro- 
dynamic simulations, assuming spherically symmetric initial 
conditions. We provide fitting formulas for the mean and peak 
accretion rates, and the period between accretion rate bursts 
as a function of the parameters we explore, including radia- 
tive efficiency, black hole mass, gas density, temperature and 
spectrum of radiation. In § 2 we introduce basic concepts and 
definitions in the problem. Numerical procedures and physi- 
cal processes included in the simulations are discussed in § 3. 
Our simulation results and the parameter study are shown in 
§ 4. In § 5 we lay out a physically motivated model that de- 
scribes the results of the simulations. Finally, a summary and 
discussion are given in § 6. 

2. BASIC DEFINITIONS 

2.1. Bondi Accretion and Eddington luminosity 

The assumption of spherical symmetry allows to trea t ac- 
cretion problems analytically. The solution (Bond ill952l) pro- 
vides the typical length scale at which gravity affects gas dy- 
namics and the typical accretion rate as a function of the black 
hole mass M&ft, ambient gas density p^ and sound speed 
c s _ 00 . The Bondi accretion rate for a black hole at rest is 

M B = / hr\BrlpooC s ,oo 
C 2 M 2 

= 4^A sPoo ^^, (1) 

where r& = GMbhCj 2 ^ is the Bondi radius, and As is the di- 
mensionless mass accretion rate, which depends on the poly- 
tropic index, 7, of the gas equation of state P = Kp 1 : 



Ar = - 



3 7 



(2) 



The value of As ranges from e 3 / 2 /4 ~ 1.12 for an isothermal 
gas (7 = 1 ) to 1/4 for an adiabatic gas (7 = 5/3). 

However, a fraction of the gravitational potential energy of 
the inflowing gas is necessarily converted into radiation or 
mechanical energy when it approaches the black hole, sig- 
nificantly affecting the accretion process. Photons emitted 
near the black hole heat and ionize nearby gas, creating a hot 
bubble which exerts pressure on the inflowing gas. Radia- 
tion pressure may also be important in reducing the rate of 
gas inflow (see MBCO09). These processes may act as self- 
regulating mechanisms limiting gas supply to the disk from 
larger scales and, thus, controlling the luminosity of the black 
hole. We quantify the reduction of the accretion rate with re- 
spect to the case without radiative feedback by defining the 
dimensionless accretion rate 



M 
Mb~' 



(3) 



where Mr is the Bondi accretion rate for isothermal gas ( 
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M B = e 3 / 2 irG 2 M 2 



bhP°° C s 



,). This definition of X ra d is con- 



sistent with the one adopted by MBCO09. 

2.2. Luminosity and Radiative efficiency 

The Eddington luminosity sets an upper limit on the lu- 
minosity of a black hole. In this limit the inward gravita- 
tional force on the gas equals the radiation pressure from pho- 
tons interacting with electrons via Compton scattering. Al- 
though this limit can be evaded in some special cases, ob- 
servations suggest that black hole and SMBH luminosity is 
sub-Eddington. The Eddington luminosity is thus, 



LEdd = ' 



4irGM, 



bhm p C 



3.3 x ItfLr. 



M, 



hi, 



(4) 



cr T ~ \100M Q/ 

The luminosity of an accreting black hole is related to the ac- 
cretion rate via the radiative efficiency 77: L = riMc 2 . From 
the Eddington luminosity, we define the Eddington gas accre- 
tion rate AlEdd = LEddC 2 , and the dimensionless accretion 
rate and luminosity as 



M 



M 



and I 



'J-Edd 



■lEdd 



(5) 



Hence, in dimensionless units, the bolometric luminosity of 
the black hole is I = 77777, where m is the accretion rate onto 
the black hole. Note, that our definition of MEdd is indepen- 
dent of the radiative efficiency 77. Therefore, if we impose 
sub-Eddington luminosity of the black hole, the dimension- 
less accretion rate ranges between < m < -. The radia- 
te — Tj 

tive efficiency, 77, depends on the geometry of the accretion 
disk and on 777. For a thin disk, r\ ~ 0.1, whereas 77 oc m 
for a n advection dominated thick disk or f or spherical accre- 
tion (lShapirolll97l iPark & Ostrikerll2001l) . In this study we 
consider two idealized cases for the radiative efficiency. The 
case of constant radiative efficiency 77 = const; and the case 
in which the radiative efficiency has a dependence on the di- 
mensionless accretion rate and luminosity: 77 = const for 
I > 0.1 and 77 oc 777 for I < 0.1. The second case we explored 
accounts for the lower radiative efficiency expected when the 
accreted gas does not settle into a thin disk. In both formula- 
tions the the radiative efficiency is one of the free parameter 
we allow to vary and we do not find important differences be- 
tween the two cases. Observations of Sgr A*, the best studied 
case of low accretion rate onto a SMBH, suggest that the ra- 
diative efficiency is indeed low but not as low as implied by 
the sc aling 77 oc 777. Recent theoretical work by (ISharma et alj 
12007ft demonstrates that there is indeed a floor on the radiative 
efficiency. 

Because the Bondi rate, Mb does not include radiation 
feedback effect, it provides an upper limit on the accretion 
rate from large scales to radii near the black hole. The Ed- 
dington rate provides the maximum accretion rate onto the 
black hole, limited by radiation feedback at small radii. Thus, 
numerical simulations are necessary to obtain realistic esti- 
mates of the accretion rates. If the accretion rate onto the 
black hole is lower than the gas accretion from large scales, 
the accreted material accumulates near the black hole, creat- 
ing a disk whose mass grows with time. We cannot simulate 
such a scenario because it is too computationally challenging 
to resolve a range of scales from the Bondi radius to the ac- 
cretion disk in the same simulations. Here we assume that 
accretion onto the black hole is not limited by physical pro- 
cesses taking place on radial distances much smaller than the 



sonic radius. For instance, even if angular momentum of ac- 
creted gas is small and the circularization radius r c « r s , fur- 
ther inflow will be slowed down with respect to the free-fall 
rate. The rate of inflow will be controlled by angular momen- 
tum loss (e.g. torques due to MHD turbulence) and there will 
be a delay between the accretion rate at the inner boundary 
of our simulation (r m j„) and the accretion luminosity associ- 
ated with it. The effect of the aforementioned time delay on 
the feedback loop is not considered in this paper but will be 
considered in future works. We also assume that the effect 
of self-gravity is negligible in our simulations since we have 
estimate that the mass within the HII region around the black 
hole is smaller than the black hole mass for Mf,h < 1000 M . 

If the rate of gas supply to the disk is given by the Bondi 
rate, accretion onto the black hole is sub-Eddington for black 
hole masses 



M bh < 



MMeT^n^rjIi., 



Gn H ,ooVTCT] 



(6) 



where we use the notations of = Too/(10 4 K), Uh,5 = 
77# iOO /(10 5 cm -3 ) and 77_i = 7//10 -1 . Thus, in this regime 
we may assume that the accretion is quasi-steady in the sense 
that the mean accretion rate onto the black hole equals the gas 
supply from large scales when the accretion rate is averaged 
over a sufficiently long time scale. 

3. NUMERICAL SIMULATIONS 
3.1. ZEUS-MP and Radiative Transfer Module 

We perform a set of hydrodynamic simulations to under- 
stand accretion onto IMBHs regulated by radiation feedback. 
Numerical simulations of radiative feedback by black holes 
are challenging because they involve resolving a large dy- 
namic al range in leng th scales. In this study we use ZEUS- 
MP (lHayes et al.l 120061) . a modified parall el version of the 
non-r elativistic hydrodynamics code ZEUS (IStone & Normanl 
|1992|). For the present work we add a radiative transfer mod- 
ule dRicotti, Gnedin, & Shullll2001l) to ZEUS-MP to simulate 
radiative transfer of UV and X-ray ionizing photons emitted 
near the black hole. A detailed description of the numerical 
methods used to solve radiative transfer and tests of the code 
are presented in the Appendix. 

As X-ray and UV photons ionize the surrounding medium, 
different reactions take place depending on the density 
and composition of the gas. Photo-ionization changes 
the ionization fraction of H and He. The detailed evo- 
lution of the Stromgren sphere depends on the cooling 
function A(T, Z) of the gas and thus on the metallicity, 
Z, and the fraction of gas in the molecular phase. For a 
gas of primordial composition, the cooling rate depends 
on the formation rates of H~ and H2, which depend on 
both the redshift and the intensity of the local dissoci- 
ating backgroun d in t he H2 Lyman-Werner bands (e.g. 
Shapiro & Kand 119871: lAbek Annin o s. No rman, & Zhang 
1998t iRicotti. Gnedin. & Shulll I2002allbh . In addition, the 
cooling function may depend on redshift due to Compton 
cooling of the electrons by CMB photons. We adopt atomic 
hydrogen cooling for temperatures T > 10 4 K, and use a 
simple parametric function to model complicated cooling 
physics of gas at T < 10 4 K. Thus, the temperature structure 
inside the ionized bubble is appropriate only for a low metal- 
licity gas. For a subset of simulations we also include the 
effect of helium photo-heating and cooling. We assume that 
gas cooling at temperatures below is negligible in order 
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to achieve thermal equilibrium in the initial conditions far 
from the black hole. For the parameter space in which we can 
neglect the effect of radiation pressure we find (see §5.1) that 
the accretion rate is a function of the temperature both outside 
and inside the HII region. The temperature outside the HII 
region depends on the cooling function of gas at T < 10 4 K 
and on the heating sources. The temperature inside the HII 
region depends on the spectrum of radiation and cooling 
mechanism of gas at T > 10 4 K. Thus, it depends on the 
gas metallicity and the redshift at which Compton cooling 
might become important. However, for the parameter space 
we have explored we find that Compton cooling has a minor 
effect on the temperature inside and outside the Stromgren 
sphere. 

The gas heating rate depends on flux and spectral energy 
distribution (SED) of the radiation emitted near the black 
hole. We assume a luminosity of the black hole I = tym 
(see § 2.2), where m is calculated at the inner boundary in 
our simulation (typically r TO j„ ~ 3 x 10~ 5 pc). We adopt a 
single power law v~ a for the SED, where the spectral index 
a is one of the parameters we vary in our set of simulations. 

We use an operator-split method to calculate the hydrody- 
namic step and the radiative transfer and chemistry steps. The 
hydrodynamic calculation is done using ZEUS -MP, then for 
the radiative transfe r calculation we use a ray tracing module 
dRicotti et al.ll2.Q0lT) . The radiative transfer module calculates 
chemistry, cooling and heating processes along rays outgoing 
from the central black hole, and thus is easily parallelized in 
the polar angle direction. 

We perform ID and 2D simulations in spherical coordi- 
nates. In both cases we use a logarithmically spaced grid in 
the radial direction typically with 256 to 512 cells to achieve 
high resolution near the black hole. The size ratio between 
consecutive grids is chosen according to the free parameters 
of the simulation to resolve the ionization front and resolve 
the region where the gas is in free fall. In the 2D simula- 
tions we use evenly spaced grids in the polar angle direc- 
tion and compute radiative transfer solutions in each direction. 
Flow-out inner boundary conditions and flow-in outer bound- 
ary conditions are used in the radial direction (r), whereas in 
polar angle directions (9), reflective boundary conditions are 
used. 

To determine the optimal box size of the simulations we 
make sure that we resolve important length scales in the prob- 
lem: the inner Bondi radius, rb^n, the outer Bondi radius, 
i*b,oo, the sonic radius, r s and the ionization front, R s . We se- 
lect the value of the inner boundary (typically ^3 x 10~ 5 pc 
for Mbh — 100 M©) to be smaller than the sonic point or the 
inner Bondi radius (both still far larger than the Schwarzschild 
radius of the black holes). We find that once the sonic radius is 
resolved, reducing the inner boundary box size does not create 
significant differences in the results. In most cases the ioniza- 
tion front is located outside of the outer Bondi radius and the 
box size is selected to be large enough to cover both length 
scales. We select a box size that achieves the highest possible 
resolution with a given number of grids, making sure that the 
physical quantities around boundaries remain constant during 
the simulations. The box is sufficiently large to minimize the 
effect of spurious wave reflections at the outer boundary. 

In this paper, the first of a series, we adopt idealized initial 
conditions of uniform density and temperature, zero velocity 
and zero angular momentum of the gas relative to the black 
hole. In future work we will relax some of these assumptions 



by adding turbulence in the initial condition and considering 
the effect of black hole motion with respect to the ambient 
medium and considering the effect of a time-delay between 
the accretion rate at the inner boundary of our simulations 
and the accretion luminosity. We assume monatomic, non- 
relativistic ideal gas with 7 = 5/3 which is initially neutral 
(electron fraction x e ~ 10~ 5 ). In this paper we also neglect 
the effect of radiation pressure. Our goal is to add to the sim- 
ulations one physical process at a time to understand which 
feedback loop is dominant in a given subset of the parame- 
ter space. We take this approach to attempt an interpretation 
of the simulation results in the context of a physically moti- 
vated analytical description of the accretion cycle. We will 
explore the effect of radiation pressure due to HI ionization 
and Lyman-alpha scattering in future works. However, a sim- 
ple inspection of the relevant equations suggests that radiation 
pressure is increasingly important for large values of the am- 
bient gas density (uh,oo ~ 10 7 cm -3 , see MBCO09 and § 6) 
since accretion rate approaches Eddington limit. 

4. RESULTS 

4.1. Qualitative description of accretion regulated by 
radiative feedback 

Our simulations show that UV and X-ray photons modify 
the thermal and dynamical structure of the gas in the vicinity 
of the Bondi radius. A hot bubble of gas is formed due to 
photo-heating by high energy photons and sharp changes of 
physical properties such as density, temperature, and ioniza- 
tion fraction occur at the ionization front. Figure [T] shows 8 
snapshots from one of our 2D simulations. Top half of each 
snapshot shows the gas density and the bottom half shows 
the hydrogen ionization fraction. We show the periodic os- 
cillation of the density and the ionization fraction from a 2D 
simulation in FigureQ] The time evolution of the density, tem- 
perature and ionization fraction profiles for the ID simulation 
are shown in Figure [2] We can identify 3 evolutionary phases 
that repeat cyclically: 

1) Once the Stromgren sphere is formed, it expands and 
the gas density inside that hot bubble decreases maintaining 
roughly pressure equilibrium across the ionization front. At 
the front, gas inflow is stopped by the hot gas and the average 
gas density inside the bubble decreases due to the following 
two physical processes. First, the black hole continues accret- 
ing hot gas within an accretion radius, r acc , defined as the ra- 
dius where the gravitational force of the black hole dominates 
the thermal energy of the hot gas. The accretion radius is 
similar to the Bondi radius defined by the temperature inside 
Stromgren sphere, but there exists a difference between them 
since the kinematic and thermal structure of gas is modified 
significantly by the photo-heating and cooling. Second, the 
gas between r acc and the ionization front moves towards the 
ionization front due to pressure gradients. The left panel of 
Figure [3] shows inflowing gas within r acc and outflowing gas 
outside r acc . A dense shell forms just outside the ionization 
front. Thus, the mass of the shell grows because gravity pulls 
distant gas into the system at the same time that gas within the 
the hot bubble is pushed outwards. 

2) As the average density inside the hot bubble decreases, 
the accretion rate diminishes. During this process the ra- 
dius of the Stromgren sphere remains approximately constant 
since the reduced number of ionizing UV and X-ray photons 
is still sufficient to ionize the rarefied hot bubble. Figure [5] 
illustrates this. Thus, the average gas temperature, ionization 
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FIG. 1 . — Evolution of the gas density and ionization fraction in a simulation of an accreting black hole of mass Af^/i = 100 Mq, gas density nrr, oc = 10 
cm -3 , and temperature Too = 10 4 K. In each panel the top halves show the density (number of hydrogen atoms per cm 3 ) and the bottom halves show the 
ionization fraction, x e = n e /riH, of the gas. The evolutionary sequences are shown in a clockwise direction. Top panels from left to right : A Stromgren sphere 
forms fueled by ionizing photons as the black hole accretes gas. The higher pressure inside the Stromgren sphere stops the gas inflow while the black hole at the 
center consumes the hot gas inside the ionization front. Inflowing gas accumulates in a dense shell outside the hot bubble while exponential decay of the accretion 
rate occurs due to decreasing density inside the hot bubble as gas depletion continues. Although the number of emitted ionizing photons decreases, the ionized 
sphere maintains its size because of the decrease in density inside the hot bubble. Bottom panels from right to left: The density of hot gas inside the Stromgren 
sphere keeps decreasing until pressure equilibrium across the front can no longer be maintained. Middle left: The dense shell in front of the Stromgren sphere 
collapses onto the black hole and this leads to a burst of accretion luminosity. Top left: The Stromgren sphere reaches its maximum size and the simulation cycle 
repeats. 



fraction and the size of the HII region remain constant. As the 
accretion rate increases during the burst, it produces a rapid 
expansion of the Stromgren sphere radius. During one cycle 
of oscillation, there are small peaks in the Stromgren sphere 
radius which are associated with minor increases in the accre- 
tion rate. Rayleigh-Taylor (RT) instabilities develop quickly 
when the accretion rate increases. In these phases, the accel- 
eration of the dense shell is directed toward the black hole, so 
the dense shell, supported by more rarefied gas, becomes RT 
unstable. 

3) As gas depletion continues, the pressure inside the hot 
bubble decreases to the point where equilibrium at the ion- 
ization front breaks down. The outward pressure exerted by 
the hot bubble becomes too weak to support the gravitational 
force exerted on the dense shell. The dense shell of gas col- 
lapses toward the black hole, increasing dramatically the ac- 
cretion rate and creating a burst of ionizing photons. The ion- 



ization front propagates outward in a spherically symmetric 
manner, creating a large Stromgren sphere and returning to 
the state where the high pressure inside the Stromgren sphere 
suppresses gas inflow from outside. 

4.2. Comparison of ID and 2D simulations 

In agreement with previous studies, our simulations show 
that radiation feedback induces regular oscillations of the ac- 
cretion rate onto IMBH. This result is in good agreement with 
numerical work by MCB09 for accretion onto a 100 M Q black 
hole from a high density (nn.oo — 10 7 cm~ 3 ) and high tem- 
perature (Too = 10 4 K) gas. Periodic oscillatory behavior is 
found in all our simulations for different combinations of pa- 
rameters, when assuming spherically symmetric initial condi- 
tions and a stationary black hole. This oscillation pattern is 
quite regular and no sign of damping is observed for at least 
~ 10 cycles. 
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10-* 10- 3 0.01 0.1 
r(pc) 

FIG. 2. — Top to bottom : Radial profiles of density, temperature and 
neutral/ionization fractions in ID simulation for »7=0.1, Af(,^=100 Mq, 
Tin oo=10 5 cm -3 and Too=10 4 K. Different lines indicates profiles at dif- 
ferent times: t=0.0 (dotted), t=1.13xf0 4 (solid), t=1.28x 10 4 (short dashed), 
t=1.43xl0 4 (long dashed), t=1.58xl0 4 (dot-short dashed), t=1.71xl0 4 yr 
(dot-long dashed). Solid lines : at the maximum expansion of the Stromgren 
sphere. Dot-long dashed lines : at the collapsing phase of dense shell. Physi- 
cal properties inside the Stromgren sphere change as a function of time. The 
number density and temperature of hydrogen decrease with time after the 
burst. The neutral fraction increases as a function of time from the burst. 

For the same parameters, our ID and 2D simulations are 
nearly identical in terms of oscillatory behavior in accretion 
rate and Stromgren sphere size. Figure|4]shows accretion rate 
in ID and 2D simulations for M bh =\QQ M Q ,T 00 = 10 4 K and 
n H,oo — 10 5 cm~ 3 . Note the similar pattern in accretion rate 
and period between bursts. This indicates that the ID result 
adequately represents 2D cases when the accretion flow does 
not have significant angular momentum. 

Moreover, this result demonstrates that RT instabilities 
which we observe in the 2D simulations do not affect the 
mean accretion rate or the period of oscillations. The RT in- 
stability develops during the phase when the dense shell in 
front of the ionization front is supported against gravitational 
accretion by the low densit y medium inside the hot bubble 
dWhalen & Normanll2008bllah . The top panels in Figure Q] 
show small instabilities when ionization fronts move outward, 
which largely decay over time. The pressure gradient inside 
the Stromgren sphere creates an outward force which helps 
suppress the development of the instability. 

In summary, we believe that ID simulations can be used 
in place of higher dimension simulations to determine the cy- 
cle and magnitude of the periodic burst of gas accretion onto 
IMBH. This allows us to reduce the computational time re- 
quired to explore a large range of parameter space. 

4.3. Parameter space exploration 

In this section we present the results of a set of ID sim- 
ulations aimed at exploring the dependence of the accretion 
rate and the period of oscillations of the black hole luminos- 
ity as a function of the black hole mass, Mbh, the ambient gas 
density, nn,oo, temperature, X^, and the radiative efficiency 
77. In § 5 we present results in which we allow the spectrum 
of ionizing radiation to vary as well. The accretion can be 
described by three main parameters: T cyc i e , the mean period 
between bursts, \ r ad,max, the maximum value of the dimen- 
sionless accretion rate (at the peak of the burst), and (X ra d), 
the time-averaged dimensionless accretion rate. These param- 



eters are typically calculated as the mean over ~ 5 oscillation 
cycles and the error bars represent the standard deviation of 
the measurements. 

After reaching the peak, the luminosity decreases nearly ex- 
ponentially on a time scale r on , that we identify as the dura- 
tion of the burst. Both r on and the duty cycle, fduty, of the 
black hole activity (i.e., the fraction of time the black hole is 
active), can be expressed as a function of T cyc i e , X ra d,max and 

(X r ad)- 



_ {Xrad) 
Ton = Tcycie 



fdu 



ty - 



^rad.rnax 
Ton {^rad) 



Tcycie 



A 



rad,max 



(7) 
(8) 



The values of A 



and fduty as a function of the black 



rad,max 

hole mass, the density and the temperature of the ambient 
medium are important for estimating the possibility of de- 
tection of IMBHs in the local universe because these values 
provide an estimate of the maximum luminosity and the num- 
ber of active sources in the local universe at any time. On 
the other hand, the mean accretion rate, (X ra d) is of critical 
importance for estimating IMBH growth rate in the early uni- 
verse. 

The four panels in Figure [7] summarize the results of a set 
of simulations in which we vary the free parameters one at a 
time. We find that, in most of the parameter space that we 
have explored, the period of the oscillations and the accre- 
tion rates are described by a single or a split power law with 
slope j3. In the following paragraphs we report the values of f3 
derived from weighted least squares fitting of the simulation 
results. The weight is 1/cr where a is the standard deviation 
of (X ra d) or Xrad,max over several oscillations. 

a) Dependence on the radiative efficiency 

First, we explore how the accretion depends on the radiative 
efficiency r\. This parameter describes the fraction, 77, of the 
accreting rest mass energy converted into radiation while the 
remaining fraction, 1 — 77, is added to the black hole mass. We 
have explored both constant values of the radiative efficiency 
and the case 77 oc m for I < 0.1 (see § 2.2). The simulation 
results shown in this section are obtained assuming 77 is con- 
stant. We find similar results for (X ra d), X ra d,max and T cvc ie 
when we assume 77 oc m. The radiative efficiency for a thin 
disk is about 10%. Here, we vary 77 in the range: 0.2% to 
10%. The other free parameters are kept constant with val- 
ues n Hoo = 10 5 cm~ 3 , M bh = 100 M and = 10 4 K. 
Figure|6]shows the accretion rate as a function of time for dif- 
ferent values of the radiative efficiency: 77 = 0.1, 0.03, 0.01 
and 0.003. Panel (a) in Figure [7] shows the dependence on 77 
of the three parameters that characterize the accretion cycle. 
The maximum accretion rate increases mildly with increas- 
ing 77 (log slope (3 — 0.13 ± 0.06). The average accretion 
rate is (X ra d) ~ 2.9% ± 0.2%, is nearly independent of 77 
(/3 = —0.04 ± 0.01). The period of the oscillations increases 
with 77 as T cyc i e oc 77 1 / 3 . We also show the simulation re- 
sults including helium photo-heating and cooling, shown as 
open symbols in the same panel of Figure [7] We find that 
including helium does not change the qualitative description 
of the results, but does offset the mean accretion rate, that is 
~ 41% lower and the period of the accretion bursts, that is 
~ 42% shorter. This offset of the accretion rate and period 
with respect to the case without helium is due to the higher 
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FIG. 3. — Gas density and velocity field for the simulation with r\ = 0.1, M b h = 100 Mq, tin oo = 10 5 cm -3 , and Too = 10 4 K. Left: When a Stromgren 
sphere is formed, gas inside the hot bubble is depleted by accretion onto the black hole and the outflow toward the dense shell due to pressure gradient. Right: 
Gas depletion inside the Stromgren sphere leads to the collapse of the dense shell, creating a burst of accretion. 
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FIG. 4. — Accretion rates as a function of time in ID and 2D simulations 
with »y = 0.1, M bh = 100 M , n H ,oo = 10 5 cm" 3 and Too = 10 4 K. 
Both results show similar oscillation patterns with the same period and aver- 
age accretion rate. 

temperature of the gas inside the HII region surrounding the 
black hole. 

b) Dependence on black hole mass 

We explore a range in black hole mass from 100 M Q to 
800 M Q , while keeping the other parameters constant (r\ 
0.1. n H:OQ = 10 5 cm" 3 and = 10 4 K). The results are 
shown in panel (b) of Figure [7] The mean accretion rate is 
(X ra d) ~ 2.7% ± 0.4% and the maximum accretion rate is 
A 



rad.max 



42% ± 12% (/3 



-0.26 ± 0.20). They are 



both independent of Mbh within the error of the fit. The pe- 
riod of the bursts is well described by a power-law relation 



T~cycle 



oc M, 



2/3 
bh 



c) Dependence on gas density of the ambient medium 



0.1 



0.01 





10 4 1.5x10" 
time (yr) 



2x10" 2.5x10" 



FIG. 5. — Evolution of Stromgren radius with time for 2D simulation with 
r\ = 0.1, M bh = 100 M , n Hj0 o = 10 5 cm" 3 and Too = 10 4 K. 
The solid line shows the mean size of the Stromgren radius and dotted lines 
show the minimum and maximum Stromgren radii. It shows the same period 
of oscillation seen in accretion rate as a function of time. In general, the 
Stromgren radius is proportional to the accretion rate which determines the 
number of ionizing photons. When the accretion rate is maximum, the size 
of the Stromgren sphere also has its maximum size. 



Panel (c) in Figure [7] shows the dependence of accretion rate 
and burst period on the ambient gas density, nn.oo- We ex- 
plore a range of n# iOC from 5 x 10 3 cm -3 to 10 7 cm -3 , 
while keeping the other parameters constant at 77 = 0.1, 
M bh = 100 M Q and = 10 4 K. For densities n H ,oo > 
10 5 cm' 3 , (Xrad) and X ra d.max are insensitive to n H ,oo W = 
-0.04 ±0.08 and fi = -0.18±0.13,respectively). However, 



rad.max 



are proportional 



for nn.oo < 10 cm 3 , (X ra d} and A 
to n 1 ^ (fi = 0.44±0.02 and /3 = 0.37±0.09, respectively). 

The bottom of Figure |7jc) shows the effect of density in 
determining the oscillation period. For densities nn.oo > 
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time(yr) 

FIG. 6. — Dependence of accretion rate and period of oscillations on the 
radiative efficiency r\. From top to bottom the evolution of accretion rate is 
shown for r] = 0.1, 0.03, 0.01 and 0.003. The peak accretion rate does 
not change much with r?, but intervals between oscillations decrease with 
decreasing r\. 

10 5 cm" 3 , T cyc i e is fitted well by a power law with T cyc i e oc 
n H<x> an< ^ f° r me densities Uh,oo < 10 5 cm -3 it is fitted 
well by a power law T cyc i e oc n^^. However, T cyc ie at 
^i/,oo = 10 7 cm -3 is lower than predicted by the power 
law fit for tih.oo > 10 5 cm -3 . Although Figure |4] [6] do not 
show clearly the magnitude of accretion rate during the inac- 
tive phase, it is evident in a log-log plot that accretion rate at 
minima is 4 orders of magnitude lower than during the peak of 
the burst. This is the case for all simulations but the ones with 
"■ff,oo = 10 7 cm~ 3 in which the accretion rate at minima is 2 
orders of magnitude higher than in all other simulations. The 
simulations show that the ambient gas density is an important 
parameter in determining the accretion luminosity and period 
between bursts of the LMBH One of the reasons is that the 
gas temperature inside the hot ionized bubble and the thick- 
ness and density of the dense shell in front of it depend on 
the density via the cooling function. The drop in the accretion 
rate we observe at low densities can be linked to an increase of 
the temperature within the sonic radius with respect to simula- 
tions with higher ambient density. This results in an increase 
in the pressure gradient within the ionized bubble that reduces 
the accretion rate significantly. 

d) Dependence on the temperature of the ambient medium 

Panel (d) in Figure [7] shows the dependence of accretion rate 
and period of the bursts on the temperature of the ambi- 
ent medium, T x . We vary T x from 3000 K to 15000 K 
while keeping the other parameters constant at r/ = 0.1, 
M bh = 100 M Q and nj^oc = 10 5 cm" 3 . We find (X ra d) and 

5/2 

Kad.max depend steeply on as TJ, (p = 2.44 ± 0.06). 
Except for the simulation with = 3000 K, the period 
of the accretion cycle is fitted well by a single power law 

rp-l/2 

Tcycie OC i oo 

5. ANALYTICAL FORMULATION OF BONDI ACCRETION WITH 
RADIATIVE FEEDBACK 



In this section we use the fitting formulas for (X ra d), 
X ra d,max and T cy cic obtained from the simulations, to for- 
mulate an analytic description of the accretion process. For 
ambient densities nn.oo > 10 5 cm~ 3 , we have found that 
the dimensionless mean accretion rate (X ra d) depends only 
on the temperature of the ambient medium. It is insensitive to 
T), and rijj )00 . Thus, for nn,oo > 10 5 cm~ 3 we find 

(X rad ) ~ 3.3% Till n H °f 4 ~ 3.3% T 5 J% (9) 

while for njj,cx> < 10 5 cm~ 3 we find 

(X rad ) ~ 3.3% I^g n#J. (10) 

As mentioned above, the dependence of (X rad ) on the den- 
sity is due to the increasing temperature inside the ionized 
bubble at low densities. The period of the accretion cycle de- 
pends on all the parameters we have investigated in our simu- 
lation. In the range of densities nn.oc > 10 cm -3 we find 

Tcycie = (6 x 10 3 yr) M* (2 n"* T~J 4 (11) 

where we use the notation of Mbh,i = Mjy,/(10 2 M©). How- 
ever, at lower densities nn,oo < 10 5 cm~ 3 , we find 

rcyde = (6 x 10 3 yr) ^ M% i 2 n"* T~$ A (12) 

in which only the dependence on tih,5 changes. The different 
dependence of T cy cie on nn.oo is associated with a change 
of the mean accretion rate (X ra d) for each density regime. 
The deviation of T cyc ie from the power law fit at n# i00 = 
10 7 cm -3 is not associated with any variation of the mean 
accretion rate. Our value of T cyc ie for uh,oo = 10 7 cm~ 3 is 
in good agreement with the value found by MCB09. 

5.1. Dimensionless accretion rate : {X rad ) 

In this section we seek a physical explanation for the re- 
lationship between the mean accretion rate {X rad } and the 
temperature of the ambient medium found in the simulations. 
The model is valid in all the parameter space we have ex- 
plored with a caveat in the low density regime {uh.oo < 3 x 
10 5 cm -3 ) and at low ambient temperatures (X^ < 3000 K). 

Figure [8] shows the time-averaged temperature profiles for 
simulations in which we vary 77, My,, Uh,oo and T^. In the 
case of different My, the radii are rescaled so that direct com- 
parisons can be made with the case of 100 M . Vertical lines 
indicate the accretion radius r acc , inside of which gas is ac- 
creted and outside of which gas is pushed out to the ionization 
front. We find that the value of r acc is generally insensitive to 
the parameters of the simulation as is the gas temperature at 

^ace- 
Accretion onto the black hole of gas inside the hot ionized 
sphere is limited by the thermal pressure of the hot gas and by 
the outflow velocity of the gas that is produced by the pressure 
gradient inside the Stromgren sphere. Thus, the accretion ra- 
dius, r acc , is analogous to the inner Bondi radius, r&,i n , mod- 
ified to take into account temperature and pressure gradient 
inside the hot bubble. 

Let us assume that the average accretion rate onto the black 
hole is 

(M) = 4:irX B rl C cPinC s ,im ( 13 ) 

where pi n and c s ^ n (and the corresponding temperature Ti„) 
are the density and the sound speed at r acc - Based on the 
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FIG. 7. — For each panel, peak accretion rate, average accretion rate and period between bursts are shown from top to bottom as a function of a given parameter. 
Error bars represent the standard deviation around the mean values over ~ 5 accretion cycles, (a) Dependence on r\. (X ra d) ~ const while r cy 
Open symbols indicate the simulations including helium photo-heating and cooling, which show ~ 41% lower accretion rate and ~ 42% shorter period, (b) 
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results illustrated in Figure [8] we expect the accretion rate to 
depend only on pi n , since r acc and c s ^ n can be taken to be 
constants. 

When a Stromgren sphere is formed, the gas inside the hot 
bubble expands and its density decreases. Inside the ioniza- 
tion front the temperature is about 10 4 — 10 5 K. Thus, assum- 
ing pressure equilibrium across the ionization front we find 
the dependence of pi n on T^: 



of gas. At low densities, r cyc i e oc 
Average accretion rate (X ra d) T^ 5 . With an exception at lowest temperature 

We find / = r acc /rb t i n ~ 1.8 and the temperature at r acc is 
Ti n ~ 4 x 10 4 K independent of rj, Mbh, n-H,oo and for a 
fixed spectral index of radiation a = 1.5. The dimensionless 
accretion rate inside of the Stromgren sphere normalized by 
the Bondi accretion rate in the ambient medium is then : 



(A r ad) — Aj 



2 

r accPin c s,in 



Tel 

' T, 

-*- 1.1 



(14) 
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FIG. 8. — Time-averaged temperature profiles as a function of simulation parameters. Different vertical lines indicate accretion radii, r a cc, for each parameter. 
Top left : r) ranges from 0.1 {solid line) to 0.003 (long dashed line). Top right : ranges from 100 Mq (solid) to 800 Mq (long dashed). Bottom left : Density 
ranges from 10 4 cm -3 (solid) to 10 7 cm -3 (long dashed). Bottom right : Too ranges from 5000K (solid) to 15000K (long dashed). 



(15) 



5.1.1. Dependence on temperature at accretion radius 



where we have used Xb = 1/4 appropriate for an adiabatic 

5 /2 

gas. Thus, (Xrad) oc Too which is in agreement with the sim- 
ulation result, given that r acc and Tj„ remain constant when 
we change the simulation parameters. However, r acc and Ti„ 
may not stay constant if we modify the cooling or heating 
function, for instance by increasing the gas metallicity or by 
changing the spectrum of radiation; this result suggests that 
the accretion rate is very sensitive to the details of the tem- 
perature structure inside the Stromgren sphere which shows a 
dependence on Hh,oo- The temperature profile changes sig- 
nificantly for riH,oo < 3 x 10 4 cm" 3 and this is probably the 
reason why our model does not fit perfectly (X ra d) from the 
simulations in the lower density regime. In the next section 
we test whether Equation dT3l l is still a good description of 
our results when we change the thermal structure inside the 
HII region. 



In this section we study the dependence of the accre- 
tion rate on the time-averaged temperature Tj„ at r acc . We 
change the temperature Tj„ by varying the spectral index a 
of the radiation spectrum and by including Compton cool- 
ing of the ionized gas by CMB photons. Here we explore 
the spectral index of the radiation spectrum in the range a — 
0.5, 1.0, 1.5, 2.0, 2.5 with the energy of photons from 10 keV 
100 keV. The other parameters are kept constant at r) = 0.1, 
M bh = 100 M©, n H ^ = 10 5 cm" 3 and = 10 4 K. 

Figure [9] shows the different time-averaged temperature 
profiles for different values of a. Spectra with lower val- 
ues of the spectral index a produce more energetic photons 
for a given bolometric luminosity, increasing the temperature 
inside the ionized bubble. Simulations show that the aver- 
aged accretion rate (X ra d) increases for softer spectrum of 
radiation. Different slopes (0.5 < a < 2.5) of the power- 
law spectrum lead to different T m (59000 K to 36000 K) and 
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FIG. 9. — Le/f : Average temperature profiles of the HII region as a function of spectral index a. Smaller a results in smaller r acc and higher Ti n . Right : 
Relation between temperature at r acc and average accretion rate (\ ra d)- We find {\ ra d) °c T~ n . 



(Kad) (0.0076 to 0.0509). Adopting a harder spectrum (with 
a = 0.5) instead of the softer (a = 2.5) increases Ti n by a 
factor of 1.6 and (X r ad) decreases by a factor 6.7. The fit to 
the simulation results in Figure [9] show that (\ ra d) depends 
on temperature at r acc as 



(Kad) oc T m 4 oc c s 



(16) 



The dependence on c s .i n differs from equation $15[ . How- 
ever this is not surprising because in these simulations the val- 
ues of r acc and c Si i n do not remain constant while we vary the 
value of the spectral index a. This is due to a change of the 
temperature and pressure gradients within the HII region. The 
accretion radius, r acc , can be expressed as a function of the 
Bondi radius inside the hot bubble, r^in = GM hC~\ n . From 
the simulations we obtain the following relationship between 
these two radii: 



We estimate f duty by comparing \ ra d,max and (A rod ) using 
Equation (HJ. This quantity gives an estimate of what fraction 
of black holes are accreting gas at the rate close to the max- 
imum. Within the fitting errors, the log slopes of X ra d,max 
and (X ra d) as a function of the parameters M h, are zero. 
Thus, we assume that the dimensionless accretion rates are 
independent of these parameters. 

For n,H,oo > 10 5 cm" 3 , \ ra d,max can be expressed as 
Kad.max ~ 0.55 rf_:\ 3 n^° 5 18 T^° 4 and the dependence of 
fduty on these parameters can be expressed using equation 
© as 



duty ' 



'6% 7? 



-0.13 n 0.14 



n0.5 



H,5 1 oo,4 



(19) 



/ = 



T 



4 x 10 4 K 



-0.7±0.2 



(17) 



where we include the mild dependence of (\ ra d) on 
the density. For n H ,oo < 10 5 cm -3 , \ ra d, m ax ~ 
0.55 T)_l n jfl T^f 4 has a different power law dependence 
on the density and we get fduty as 

f „-0.13 „,0.07 

J duty ~D/o T]_ 1 



rpO.5 

H,5 1 ooA 



(20) 



Thus, if our model for the accretion rate summarized by 
equation (flj) is valid, we should have: 



(Kad) — ~. 



1 r l 



accPin c s,in 
4 r b,oo^ 



(l-8) s 



C - 5 ' 9 c 3 

s,tn s,oo 



T 



4 x 10 4 K 



(18) 



in agreement with the simulation results (\ rad ) oc T^T^ 4 
where the dependence on T- m was not explored initially. Thus, 
Bondi-like accretion on the scale of r acc is indeed a good ex- 
planation of our results. Given the steep dependence of the 
value of accretion rate (X ra d) on Ti n it is clear that it is very 
sensitive on the details of the thermal structure inside the HII 
region. This means that (X r ad) depends on the spectrum of 
radiation and gas metallicity. 



5.2. Accretion rate at peaks and duty cycle: A 



rad.max- 



f 



duty 



where fduty shows a milder dependence on the gas density. 
Thus, we expect about 6% of IMBHs to be accreting near the 
maximum rate at any given time. This value depends weakly 
on r\, n H .oa and T^. 

5.3. Average period between bursts : T cyc i e 

In this section we derive an analytical expression for the 
period of the luminosity bursts as a function of all the pa- 
rameters we tested. Although T cyc i e shows a seemingly com- 
plicated power law dependencies on the free parameters, we 
find that T cyc ; e is proportional to the time-averaged size of the 
Stromgren sphere. This is shown in Figure [TU] The linear re- 
lation between T cyc i e and the average Stromgren radius (R s ) 
explains the dependence of r cyc i e on every parameter consid- 
ered in this work. 

The number of ionizing photons created by accretion onto 
a black hole is determined by the average accretion rate and 
the radiative efficiency ?/. The average accretion rate itself can 
be expressed as a fraction of the Bondi accretion rate (X rad )- 
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Therefore, the average number of ionizing photons emitted 
near the black hole can be expressed as 



N ion ocrj{\ rad )M B 

oc rj{\ rad ) — 5 Poo • 



(21) 



It follows that: 



T 'cycle ' 



(Rs) 



oc 



Vout 
1 

3 



3iVi, 



47rav ec ra^ 



ri(Xrad) , 



/Ooo 



(22) 



where we find u ou 



3 ^s.zn ■ 



Ignoring constant coefficients 



and using equation (0 for nu oo > 10 5 cm 3 , we find : 

i i—i —i 

Tcycle OC T] s M b \n 2 , (23) 



or using equation ( TTOb for riff l0O < 10 5 cm 3 , we find : 

1 2^ -L _ 1 

Tcycle OCT) 3 M^n^T^ 



(24) 



which are exactly as in the empirical fitting formulas in both 
density regimes and also in good agreement with the ana- 
lytical work by MBCO09. This explains the dependence of 
T cyc i e on any tested parameter t], Mbh, n>H,oo an d ?oo- In Fig- 
ure[lO]we also show simulation results assuming r\ oc m. All 
simulations show the same relationship between T cyc i e and 
(R s )- However, the simulation with the highest ambient den- 
sity (jih,oo = 10 7 cm -3 ) deviates from the linear relation- 
ship, but is in agreement with the numerical simulation by 
MCB09. It appears that in the high density regime T cyc ie de- 
creases steeply with decreasing (R s ). 

We can interpret T C y C ie as the time scale at which the gas 
inside HII region gets depleted. If the gas depletion inside the 
Stromgren sphere is dominated by the outward gas flow, then 
Tcycie oc (R s ) / ' c s i n in agreement with the empirical linear 
relation in Figure [10] However, the depletion time scale may 
be different if the accretion by the black hole dominates gas 
consumption inside the Stromgren sphere. We can derive this 
time scale as 
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Roughly, we expect T cyc ie = mm(t out ,t m ). So, for 
(i? s )/r occ < 3, the period of the cycle scales as (R s ) 3 - This 
may explain the deviation of the period for njjoo = 10 7 cm -3 
from the linear relation. We see in Figure [8] that the ratio 
(i? s )/r occ ~ 5 for nn,oo — 10 7 cm -3 which is much smaller 
than the ratio found for other densities. 

5.3.1. Rayleigh-Taylor instability 

In 2D simulations we find that RT instability develops 
across the Stromgren radius, but it decays on short time scales. 
This can be explained by the pressure gradient inside the 
Stromgren sphere which does not allow the RT grow. In the 
linear regime the growth time scale of the RT instability of 
wavelength A is 
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FIG. 10. — Period of accretion bursts as a function of the average Stromgren 
radius. All simulation results from all the parameters are plotted together. 
The average size of the Stromgren sphere shows a linear relation with the 
period T cyc i e . The only exception happens at the highest density (nj iCO = 
10 7 cm -3 ), but this result is in agreement with the work by MCB09 (symbol 
with error bar) 

where p s h is the density of the shell and g ~ GMbh(Rs)~ 2 
is the gravitational acceleration at the shell radius. Thus, RT 
timescale can be expressed as 
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So during one cycle perturbations grow on scales: 
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where r^.in is me inner Bondi radius. Thus only instability on 
angular scales 9 ~ \ rt /2k(R s ) < r btin / (2-k) 2 (R s ) grow in 
our simulation. 

6. SUMMARY AND DISCUSSION 

In this paper we simulate accretion onto IMBHs regulated 
by radiative feedback assuming spherical symmetric initial 
conditions. We study accretion rates and feedback loop pe- 
riods while varying radiative efficiency, mass of black hole, 
density and temperature of the medium, and spectrum of radi- 
ation. The aim of this work is to simulate feedback-regulated 
accretion in a wide range of the parameter space to formu- 
late an analytical description of processes that dominate the 
self-regulation mechanism. Thus, in this first paper we keep 
the physics as simple as possible, neglecting the effect of an- 
gular momentum of the gas, radiation pressure and assum- 
ing a gas of primordial composition (i.e. metal and dust 
free). We will relax some of these assumptions in future 
works. However, the parametric formulas for the accretion 
presented in this paper should provide a realistic description 
of quasi-spherical accretion onto IMBH for ambient gas den- 
sities tih.oo ^ 10 5 — 10 6 cm -3 , as radiation pressure should 
be minor for these densities. 

We find an oscillatory behavior of the accretion rate that can 
be explained by the effect of UV and X-ray photo-heating. 
The ionizing photons produced by the black hole near the 



Accretion onto IMBHs Regulated by Radiative Feedback 



13 



gravitational radius increase gas pressure around the black 
hole. This pressure prevents the surrounding gas from being 
accreted. An over-dense shell starts to form just outside the 
Stromgren sphere. Due to the decreased accretion rate, the 
number of emitted ionizing photons decreases and the density 
inside the Stromgren sphere also decreases with time. Gas 
accretion onto the black hole is dominant in decreasing the 
density inside the HII region only for ambient gas density 
n<H,oo ^ 10 7 cm~ 3 ; for lower values of the ambient gas den- 
sity the gas inside the HII region is pushed outward toward 
the dense shell by a pressure gradient that develops behind the 
ionization front. Eventually, the pressure gradient inside the 
Stromgren sphere is not able to support the weight of the over- 
dense shell that starts to fall toward the black hole. The accre- 
tion rate rapidly increases and the Stromgren sphere starts to 
expand again. 

However, the introduction of a small, non-zero, angular mo- 
mentum in the flow could change the time-dependent behav- 
ior of accretion and feedback loop. The inflow rate in the 
accretion disk that will necessarily develop, and that is not re- 
solved in our simulations, is typically much slower than the 
free-fall rate since the viscous time scale in units of the free- 
fall time is t V i sc /tff ~ a " 1 M 2 where a is the dimen sionless 
parameter for a thin disk (Shakura & Sunyaev 1973) and M. 
is the gas Mach number. Therefore, angular momentum may 
produce a long delay between changes in the accretion rate 
at the inner boundary of our simulation and their mirror in 
terms of output luminosity. Hence a ~ 0.01 — 0.1 and M. 
at the inner boundary of our simulations is of order of unity, 
time delays of 10-100 free-fall times are shorter if the disk is 
smaller than the inner boundary of the simulation. We have 
started investigating the effect of such a delay on the periodic 
oscillations and preliminary results show that the period of the 
oscillations can be modified by the time delay but the oscil- 
latory behavior is still present (at least for delays of 10-100 
free fall times calculated at the inner boundary). As long as 
the time delay is shorter than the oscillation period, that de- 
pends mainly on the gas density, it does not seem to affect the 
results. We are carefully investigating this in the low and high 
density regimes where the oscillation pattern and the periods 
are different. At low densities a time delay of a few hundred 
free-fall times is much smaller compared to the oscillation pe- 
riod, whereas at the high densities the maximum time delay 
that we have tested is comparable to the oscillation period. 
We will publish more extensive results on this effect in our 
next paper in this series. 

We find that the average accretion rate is sensitive to the 
temperature of the ambient medium and to the temperature 
profile inside the ionized bubble, and so depends on the gas 
cooling function and spectral energy distribution of the radi- 
ation. The period of the accretion bursts is insensitive to the 
temperature structure of the HII region, but is proportional to 
its radius. 

Our simulations show that ID results adequately re- 
produce 2D results in which instabilities often develop. 
We find that the accretion rate is expressed as X ra d — 

3% T^ 5 4 (Tin/4 x 10 4 K)~ 4 . We also derive T cycie as a 
function of r], M&fc, nn.oo and Too. The dependencies 
of (Kad) and T cyc i e on our free parameters can be ex- 
plained analytically. Assuming pressure equilibrium across 
the Stromgren sphere is a key ingredient to derive the depen- 
dence of (X ra d) on Too, whereas the linear relation between 
the average size of the Stromgren sphere and T cyc i e is used 



to derive the dependence of T cyc ie on all the parameters we 
varied. 

The qualitative picture of the feedback loop agre es with the 
description of X-ray bursters in lCowie et al.l(fl978|) . After ex- 
trapolating our analyt ical formulas to bla ck holes of a few so- 
lar masses studied by ICowie et al.1 (11978). we find that the av- 
erage accretion rate is in good agreement (L ~ 2 x 10 35 erg/s). 
However, the details of the accretion rate as a function of time, 
the burst pe riod and peak accr etion rates show qualitative dif- 
ferences. ICowie etaLl (119781) simulations do not show peri- 
odic oscillation while our simulations have well-defined fast 
rise and exponential decay of accretion followed by quiescent 
phases of the accretion rate. This regular pattern of accre- 
tion bursts is possible only when spherical symmetry is main- 
tained on relatively large scal es during o s cillations. An ax- 
isymmetric radiation s ource (Proaa 20071 iProea et all 120081: 
iKurosawa et alJ 120091: [Kurosawa & Proga 2009a.tJ or inho- 
mogeneous initial condition on scale of the Bondi radius can 
break the symmetry. 

Our simulations are also in excellent qualitative agree- 
ment with simulations by MCB09 that studied accretion onto 
100 Mq black hole for the case nu,oo — 10 7 cm" 3 . How- 
ever, we find a dimensionless accretion rate (X ra d) ~ 3% 
((Xrad) ~ 2% including helium heating/cooling) that is about 
one order of magnitude larger than in MCB09. The cycle pe- 
riod, T cyc ie, is in better agreement since T cyc ie °c (R s ) oc 
(Xrad) 1 ^ 3 - The discrepancy in the mean accretion is likely 
produced by the effect radiation pressure on HI, that becomes 
important for large n# j00 and that we have neglected. In 
addition, our results indicate that the qualitative description 
of the feedback loop starts to change at ambient densities 
> 10 6 — 10 7 cm -3 : the oscillation period decreases much 
more rapidly with increasing ambient density as gas deple- 
tion inside the ionized bubble becomes dominated by accre- 
tion onto the IMBH. The accretion luminosity during the qui- 
escent phase of the accretion also increases and the two phases 
of growth and collapse of the dense shell become blended into 
a smoother modulation of the accretion rate. Hence, further 
numerical studies are required to characterize accretion onto 
IMBH in the high-density regime. 

As mentioned above, in this study we have neglected three 
important physical processes that may further reduce the ac- 
cretion rate: 1) Compton heating, 2) radiation pressure and 
3) Lyman-a scattering processes. The importa nce of these 
proces ses is thoroughly discussed in MBCO09. iRicotti et al.l 
(2008) and MBCO09 find that Compton heating is not an im- 
portant feedback mechanism in regulating the accretion rate 
onto IMBHs in which the gas density inside the Stromgren 
sphere is roughly independent of the radius. However, 
MBCO09 suggest that both radiation pressure on HI and 
Lyman-a radiation pressure can contribute to reducing the 
accretion rate onto IMBH. At higher densities accretion be- 
comes Eddington limited ((X ra d)MB ^ MEdd)- By assum- 
ing (X r ad) ~ 1% T^ 5 4 which is suggested by simulations 
including helium and manipulating Equation (O we obtain 
Mbh,2 Too 4 Uh,5 V-i > 40 as a criteria for Eddington lim- 
ited condition. This expression predicts that radiation pres- 
sure becomes important at gas density nn,5 ~ 10 — 100 
with other parameters fixed to unity. It also implies that this 
critical density depends on the black hole mass, gas temper- 
ature and radiative efficiency. We expect radiation pressure 
to play a minor role at low densities also because the accre- 
tion luminosity is negligible during the quiescent phases of 
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accretion (between accretion bursts). However, at densities 
nu,oo > 10 7 cm~ 3 the accretion rate during the quiescent 
phases is not negligible, thus radiation pressure can be im- 
portant in this regime. The main effect of Thomson radia- 
tion pressure is to prevent the accretion luminosity to exceed 
the Eddington limit. The continuum radiation pressure due 
to HI ionization can instead be important for sub-Eddington 
luminosities, but will be strong only at the location of the ion- 
ization front during the peak of the accretion burst. We will 
present more extensive results in the next paper of this series. 

The results of this study provide a first step to estimate the 
maximum X-ray luminosity and period of oscillations of an 
accreting IMBH from a medium with given physical condi- 
tions. Hence, they may be useful for modeling detection prob- 
ability of ULX originating from accreting IMBH in the local 
universe. From the average growth rate of IMBHs accreting 
in this manner it is also possible to estimate the maximum 
masses of quasars at a given redshift starting from seed pri- 
mordial black holes. One of the main motivations of this study 
is to derive simple analytical prescriptions to incorporate 



growth of seed black holes from Population III stars into large 
scale cosmological simulation. However, before being able to 
do use these results in cosmological simulation we need to 
understand the effects of relaxing our assumption of spheri- 
cal symmetric accreti on: we need to simulate a c cretion onto 
moving black holes dHovle & Lvttletonl Il939t [Shima et al. 
119851; iRuffert & Arnettlll994tlRuffertll 19961) and use more re- 
alistic initial conditions, including gas with non-zero angular 
momentum or a multi -phase turbulent ISM (Krumholz et al. 
I2005ll2006h . 

The simulations presented in this paper were carried out 
using high performance computing clusters administered by 
the Center for Theory and Computation of the Department 
of Astronomy at the University of Maryland ("yorp"), and 
the Office of Information Technology at the University of 
Maryland ("deepthought"). This research was supported by 
NASA grants NNX07AH10G and NNX10AH10G. The au- 
thors thank the anonymous referee for constructive comments 
and feedback. 
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APPENDIX 
BASIC TESTS OF THE CODE 

We test the Bondi accretion formula using ZEUS-MP for the adiabatic indexes 7 = 1.2, 1.4 and 1.6. For a given equation of 
state, the sonic point where the gas inflow becomes supersonic must be resolved not to overestimate the accretion rate Xb- The 
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left panel of Figure QT| shows the steady accretion rate as a function of the radius at the inner boundary normalized by the Bondi 
radius. Different lines show results for 7 =1.2, 1.4 and 1.6. 
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FIG . 11. — Left : Simulated Bondi accretion rate( A g ) as a function of minimum radius with given adiabatic index 7 without radiative feedback. Dashed lines 
are analytically estimated values for each 7 =1.2, 1 .4 and 1.6. In order not to overestimate accretion rate sonic point should be resolved where the velocity of 
the inflowing gas becomes supersonic. Right : Test of Stromgren radius with given number of ionizing photons. Solid line is the prediction for the given number 
of ionizing photons from 10 40 to 10 50 s — 1 . Triangle symbols represent location where ionization fraction of njj (^hi) is 0.50. Squares are for xjji = 0.90 and 
circles are for xm = 0.99. 

We also test whether our radiative transfer module produces radii of the Stromgren spheres in agreement with the analytical 
prediction: (4-7r/3)i?;? = Ni on , where R s is the Stromgren radius and Ni on is the number of ionizing photons emitted 

per unit time. The right panel of Figure [TT] shows the test of the ID radiative transfer module without hydrodynamics. Different 
symbols indicates the radii for the different ionization fractions: x e — 0.99 (circle), 0.90 (square), 0.50 (triangle). 

RADIATIVE TRANSFER MODULE AND TIME STEPPING 

Our hydrodynamic calculation is performed using ZEUS-MP, returning the density and gas energy at each time step to the 
radiative transfer module. The op erator-splitting method is applied to mediate between hydrodynamics and radiative transfer 
with a photon-conserving method ( Wh alen & Normanll2006l) . For each line of sight radiative transfer equations are solved in the 
following order: 

1 . At the inner boundary, the average inflow mass flux M is calculated. 

2. The mass flux is converted into accretion luminosity L, and thus into the number of ionizing photons for a given radiative 
efficiency r/. 

3. The photon spectrum is determined using a power law spectral energy distribution with the spectral index a. We use up to 
300 logarithmically spaced frequency bins for photons between 10 eV up to 100 keV. 

4. The ordinary differential equation for time-dependent radiative transfer cooling/heating and chemistry of the gas are solved 
using a Runge-Kutta or Semi-Implicit solver for each line of sight with a maximum of 10% error. Photo-heating, cooling 
for a given cooling function and Compton cooling are calculated. 

5. The energy density and the abundances of neutral and ionized hydrogen are updated. 

Parallelization is easily implemented in polar angle direction because radiative transfer calculations along each ray are indepen- 
dent of one another. 

RESOLUTION STUDIES 

We perform a resolution study to confirm that the number of grid zones does not affect the results. Number of zones from 384 
to 768 are tested and they all show the similar outputs in terms of accretion rate at peaks, average accretion rate, decaying shape 
and the period between peaks . Figure Q~2] shows that the details of the accretion rate history from simulations are not identical 
but the physical quantities which we are interested in (average accretion rate, peak accretion rate and period of the bursts) do 
not show significant deviation from each other. In general, a Courant number of 0.5 is used for most simulations, but we try a 
Courant number which is one order of magnitude smaller to investigate how the results are affected by reducing the hydro-time 
step by an order of magnitude. The chemical/cooling time steps are calculated independently by the radiation transfer module. 
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FIG. 12. — Comparisons between simulations of r\ = 0.1, Mbh — 100 Mq, ?iif ; oo — 10 cm 

" 3 and Too = 10 4 K with various resolution. Solid : 384 grid 
run. Dotted : 512 grid run. Long dashed 768 grid run. Short dashed : 512 grid with Courant number of 0.05. 



